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We present a method of organizing incoherent motion of a colloidal suspension to produce 
synchronized, coherent motion. This method exploits general features of rotational response 
to time-dependent forcing, and it does not require interaction between the particles. We 
report two methods of achieving orientational alignment of an ensemble of identical col- 
loids by means of a time-dependent, but spatially uniform forcing: a) a piecewise constant 
force alternating between two directions and b) a force uniformly rotating about an axis. 
The physical origin of the forcing may be e.g., sedimentation or electrophoresis. We will 
demonstrate that these forcing methods achieve alignment both by analyzing the equations 
of motion and by simulation. We find the conditions guaranteeing alignment, discuss the 
limitations of these methods, and suggest possible applications. Examples of such forcing 
include electrophoresis and sedimentation. 

PACS numbers: 05.45.-a, 82.70.Dd, 87.50.ch 
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I. INTRODUCTION 

Phase coherence plays an major role in many 
applications. For example, in nuclear magnetic 
resonance, orientational alignment of the spins 
produces a large magnetic signal which is used as 
a characterization tool. It is essential for lasers, 
where the coherence of the emitted photons im- 
bues the light with remarkable properties. In 
this paper, we will describe a method of inducing 
coherence in the context of colloidal sedimenta- 
tion. Specifically, we will describe a method of 
orientationally aligning an ensemble of colloids 
by uniformly applying a time-dependent force. 
This alignment is desirable since once the en- 
semble is aligned, any response to subsequent 
forcing will automatically be coherent. Thus it 
is possible to manipulate the orientation of ob- 
jects in the ensemble in a controlled way. 

Our interest in this example of coherence has 
grown out of a body of recent work on the sub- 
ject of sedimentation of colloidal objects [1—1]. 
These references discuss the following situation: 
an object with arbitrary shape is immersed in 
an infinite volume of fluid, which is sufficiently 
viscous to make the Reynolds number small. Hy- 
drodynamic as well as external forces (e. g., 
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gravity) are assumed to act on the object, and 
the resulting positional and orientational evolu- 
tion of the object is determined. An object often 
exhibits a simple response to constant forcing: a 
direction in the object aligns with the force, and 
the object subsequently rotates about this axis 
with a constant angular velocity [2]. We call an 
object with this behavior "axially aligning." If 
an entire ensemble of such objects were subject 
to a force, they would all become aligned with a 
common axis. Still, the objects' orientations are 
not identical. They differ from one another by a 
rotation about the force axis. 

The results [1-1] have paved the way for fur- 
ther research. For example, when the forcing is 
electrophoretic, [5] explains how to generate two 
novel responses to the applied electric field: ro- 
tation without translation and translation per- 
pendicular to the applied field. Evidently, the 
direction of rotation depends on the handedness 
of the object. This property enables one to sep- 
arate enantiomers in a shear flow [6-11]. How- 
ever, the problem of removing the orientational 
indeterminacy present after constant forcing is, 
to our knowledge, unaddressed. 

In this paper our goal is to show how to ef- 
fect orientational alignment in an ensemble of 
colloids using time-dependent forcing. We will 



study two different examples of forcing, and 
show, by analysis and simulation, that these lead 
to alignment under the suitable hypotheses. In 
Section II, we will build the mathematical frame- 
work required to study the problem. In Section 
III we will review the response to constant forc- 
ing. In Section IV, we will apply our formalism 
to treat two forcing programs and show that they 
both have the potential to align. In Section V, 
we will present results of numerical simulations 
of the two forcing programs. In Section VI, we 
will discuss the assumptions necessary to obtain 
our results, determine the minimum object size 
necessary to suppress thermal rotational diffu- 
sion, and briefly suggest possible applications. 



II. MATHEMATICAL FORMALISM 
A. Rotational Response 

We may use the established linear-response 
formalism to describe the orientation effects of 
interest [12]. We are interested in the motion 
of an ensemble of identical objects which are 
dispersed in a fluid and subject to an external 
time-dependent force Fe(t). For definiteness, we 
assume the external force is gravity; we will dis- 
cuss generalizations to other types of forcing at 
the end of this section. We consider the regime of 
"creeping flow" [12], in which inertial forces are 
negligible and the force transmitted to a mov- 
ing object by the medium is proportional to the 
object's velocity. For a rotating rigid body, the 
transmitted force F^ and torque f^ both depend 
linearly on the center of mass velocity v and an- 
gular velocity uj. 

How do controllable forces influence an ob- 
jects orientation? Since inertial forces are negli- 
gible, any external force must be exactly opposed 
by a hydrodynamic force. For an object to gen- 
erate this hydrodynamic force, it must have the 
appropriate velocity and angular velocity. The 
orientation may be altered by prescribing the an- 
gular velocity. Thus the external force can be 
used to control the orientation. Below we look 
at each step in more detail. 

As noted above, the hydrodynamic force 
must satisfy Fh = —F^. Since rotational inertia 



is also negligible, the hydrodynamic and exter- 
nal torques must be related similarly: t^ = —^e- 
Thus the hydrodynamic forces and torques are 
easily controlled. Because these relationships are 
so simple, it is convenient to refer only to the hy- 
drodynamic force and torque, and to drop their 
subscripts, referring to them simply as F and r 
respectively. 

The connection between the hydrodynamic 
force and the object's motion is given by its 
linear-response properties. To discuss these 
properties, it is convenient to introduce some 
new notation. Let i? be a characteristic size of 
the object under consideration; for definiteness 
we take the stokes radius [12]. Following [1], we 
define two six-component vectors: the general- 
ized hydrodynamic force vector F = {F, t/R)'^ 
specifying the force and torque, and the general- 
ized velocity vector V = {V,loR)'^ specifying the 
velocity and angular velocity. We noted that it 
is possible to impose a desired F on the system. 
The object is obliged to assume the V which gen- 
erates this force. In the creeping flow regime, an 
object moving with generalized velocity V gener- 
ates a generalized hydrodynamic force given by 



F = 6Tn]RFV, 



(1) 



where P is a dimensionless symmetric 6x6 ma- 
trix. To simplify (1), we choose units where 
R = Gttt] = 1. Since we seek the motional re- 
sponse required to generate our imposed force, 
we will consider the inverse of P, denoted M [12]; 
the required V is 



V 



IF. 



(2) 



To separate the roles of force and torque, we 
distinguish the four 3x3 submatrices of M that 
give the response to v and a; to F and r by writ- 
ing M the following way: 



A T^ 

T S 



(3) 



In sedimentation, the external force may be re- 
garded as acting at a single point: the center 
of buoyant mass. We choose this forcing point 
as the origin, so that the external torque on the 
body is zero. Then the rotational motion of the 



body is completely determined by the 3x3 "twist 
matrix" T: 



uJ = TF. 



(4) 



This equation captures how an applied force 
causes the orientation to change; it will be fun- 
damental in what follows. The complexity of 
this equation lies in the fact that the tensor T, 
being a property of the object, changes when 
the object is rotated. Thus, even if the force 
is constant, the angular velocity causes T to 
change, which in turn causes the angular veloc- 
ity to change. This interplay can lead to compli- 
cated dynamics as discussed in Section III. 

The T matrix of an object depends on the 
position of the origin, and hence on the forcing 
point; this point may be varied without changing 
the object's shape e.g., by changing the mass 
distribution within it. A simple transformation 
law determines how T changes under a change 
of forcing point [1]. One location in the object 
called the "center of twist" makes T symmetric. 
Moving the forcing point sufficiently far from the 
center of twist always causes the object to be 
axially aligning [1]. 

We may extend the results of this section to 
forcings of non-gravitational origin. One exam- 
ple is to increase buoyancy effects via centrifuga- 
tion. The only difficulty with using a centrifugal 
force is its non-uniformity. Otherwise the cen- 
trifugal force is indistinguishable from a gravi- 
tational force of the same strength, so the ar- 
guments of this section are valid without mod- 
ification. Another forcing to consider is elec- 
trophoresis. In this case, the electric field's effect 
on dissolved counter-ions cause the fluid to flow. 
This flow may lead to exotic motion, such as 
translation purely perpendicular to the applied 
electric field or rotation without translation [13]. 
Nonetheless, it has been found that even in this 
case, the angular velocity of an object depends 
linearly on the applied field. Thus, provided that 
T is replaced by an effective twist matrix and F 
is replaced by E, (4) continues to hold. 

Having made these observations, we will no 
longer use any property of the object besides T, 
and we use no information about the forcing ex- 
cept the vector F (with the understanding that 
it may represent an electric field or even some 



other vector parameter). We will simply study 
the differential equation (4) for general T and 
F. We have reduced the problem to pure math- 
ematics. 



B. Equation of Motion 

In this section we will derive the equation 
of motion for T implied by (4). Equation (4) 
states that when a force is applied, the object 
will rotate with the angular velocity determined 
by T and F. Between the initial time t = 
and a time t, the object will have undergone a 
finite rotation, given by an orthogonal transfor- 
mation matrix ]R(t). When the object does ro- 
tate, the initial twist matrix Tq must transform 
into the twist matrix T at time t according to 
T = MTqIR-^. From this equation, it is clear that 
the time derivative of T can be expressed terms 
of the time derivative of M, which we now exam- 
ine. To write this derivative, we introduce new 
notation: for any vector v, we define the "cross 
product matrix" v^ to be the one that, acting 
on a vector u, produces the vector v x u. The 
matrix entries of v^ are given by 



[^ 



I ik 



^ijk'^j ■ 



(5) 



Then the time derivative of M is given by 
cD^M. From this expression we find 

t = — MTqM^ 
dt 



^TqM^ + WFqR^ 



-IXttdnT 



(6) 



= [(TF)^T], 
where to get the last line we have used (4). 

III. CONSTANT FORCING 



In this section, we review the motion result- 
ing from a constant force. In so doing, we will 
find conditions on the twist matrix which are 
necessary and sufficient to ensure that the object 



is axially aligning. This section lays the ground- 
work for the following discussion of non-constant 
forcing. 

The goal, then, is to solve the differential 
equation (6) assuming constant F. The char- 
acter of these solutions depends on the matrix 
T. There are only three classes of T, each ex- 
hibiting different behavior, of which one is axial 
alignment. 

One class consists of symmetric T matrices. 
In this case there is an analogy with classical me- 
chanics. If we take (4) and replace F and T by 
the angular momentum vector L and the inverse 
of the inertia tensor I~^ respectively, we are left 
with the equation for free rotation of a body. 
Let us recall the solution to this equation. The 
standard approach is to view the motion in the 
body frame [14]. We can describe the matrix I in 
terms of its eigenvalues Ai , A2 , and A3 , and their 
corresponding eigenvectors vi, V2, and V3. Since 
we are in the body frame, both the eigenvalues 
and the eigenvectors are constant. Let us con- 
sider the generic case where the three eigenvalues 
are distinct, and let us label the eigenvalues so 
that Ai > A2 > A3. Now the vector L is allowed 
to change with time, but its length must be con- 
stant, so that L is constrained to lie on a sphere. 
Since I is symmetric, the energy E = ^LI~^L is 
also constant in time. Thus L{t) is confined to 
lines of constant E on the sphere. If L is close to 
the vi or v^ direction, then the energy E is nearly 
extremal and L is confined to a small neighbor- 
hood. Thus we see the fixed points ±^1 and ±^3 
are (neutrally) stable. However, if L is initially 
near ±7/2, it will in general move far from this 
vector and so the fixed point ±^2 is unstable. 
Translating this back into the language of sed- 
imentation, we would say that the eigenvectors 
vi and V3 of T are neutrally stable fixed points 
for the motion of F while V2 is an unstable fixed 
point. 

Now unlike I, T need not be symmetric. If 
the antisymmetric part of T is sufficiently small, 
then the eigenvalues remain real, and one of zizvi 
and one of zizv^ become unstable [2]. For this 
second class of T, two objects will not, in general, 
align with the same axis [2]. Thus we call an 
object in this class "non- aligning" . 

Finally, a T may have an antisymmetric part 



that is too big for the second class. This T be- 
longs to the third class. In this class, T has 
only one real eigenvalue, A3, the other two being 
complex conjugates, as is the case for a com- 
pletely antisymmetric T. For this third class 
of T, the motion leads to axial alignment [2]. 
In the body frame, this means that for almost 
any initial condition, F aligns along some di- 
rection. Moreover, this direction is in the one 
dimensional A3 eigenspace. This allows us to 
fix the sign of the eigenvector v^ so that it is 
the unit vector along the direction of alignment. 
Then the only F's which do not align with V3 
are those parallel with —v^. Since the criterion 
that T has only one real eigenvalue is necessary 
and sufficient for the object be axially aligning, 
we will refer to these twist matrices themselves 
as axially aligning. The motion of F in the body 
frame for the cases of symmetric, non-aligning, 
and axially aligning twist matrices is shown in 
Fig. L 

Let us now view the rotation of an axially 
aligning T in the lab frame where F is con- 
stant. The above shows that the T matrix 
must rotate so that its V3 eigendirection becomes 
aligned with F, but the orientation of the objects 
about F are unconstrained. Therefore, if an en- 
tire ensemble of objects were subject to a con- 
stant force, the ensemble would become axially 
aligned, with each object's vs aligned with F. 
We contrast this state of axial alignment with 
complete alignment, where every object in the 
ensemble has an identical orientation. The dis- 
tinction between complete and axial alignment 
is illustrated in Fig. 2. It would be convenient 
if complete alignment could be attained by use 
of a constant force. However, since each object's 
^3 is aligned with F, (4) states that each object 
has the same angular velocity: A3-P. Thus no 
change in relative orientation among objects is 
possible, and a time-dependent forcing program 
is necessary to achieve complete alignment. 



IV. COMPLETE ALIGNMENT 

We are now ready to discuss the task 
at hand: completely aligning an ensemble of 
non- interacting identical objects using a time- 
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(b) 




FIG. 1. Illustrations of possible motions of F in the object frame for differing initial conditions and twist 
matrices. The small red, blue, and green balls show the initial value of F, and the curves give the evolution 
of F. The twist matrices used to generate (a), (b), (c) were symmetric, non-aligning, and axially aligning, 
respectively. In (a) the various motions of F lie on closed curves, analogously to the angular momentum in 
the case of a freely rotating body, (b) shows the effect of adding a small antisymmetric part to the T of (a). 
There are two stable fixed points, as shown by the red and green trajectories. The blue trajectory does not 
converge to a fixed point but instead converges to a limit cycle, (c) shows the effect of further increasing 
the antisymmetric part in (b) so that complex eigenvalues appear. Here the three point converge to a single 
fixed point. 
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FIG. 2. Two snapshots of a set of objects with differing degrees of alignment, (a) Axial alignment. There 
is an axis in the object, indicated by the black line, pointing along the same direction in each instance. 
This implies that the orientation of the cars differs only by a rotation about this axis, (b) Each car in the 
ensemble is in precisely the same orientation. 



dependent force. That is, given a twist ma- 
trix T, we wish to find a time-dependent forcing 
program F{t) such that the difference between 
almost any two orientations goes to zero with 
time. In the next paragraphs we make assump- 
tions and general remarks, then we consider two 
choices of F{t). 

First, the magnitude of F plays no essential 
role in alignment: it is evident from (4) that the 
only effect of doubling F would be to speed up 
the object's motion by a factor of two. More 
generally, applying a force of varying magnitude 
is equivalent to applying a force of unit magni- 
tude and reparameterizing time. Since this repa- 
rameterization has no bearing on whether the 
final state is aligned, we may restrict our atten- 
tion to forcings of unit magnitude. Further, we 



non-dimensionalize the force by setting this unit 
equal to one. We emphasize that the force has 
unit length by writing F in place of F. 

Next we make the assumption that T is ax- 
ially aligning. We saw in the previous section 
that the motion is simpler in this case, and all 
that needs to be done is to remove a rotational 
phase indeterminacy. On the other hand, if 
the object were non-aligning, the response to a 
constant force would be more complicated: As 
shown in Fig. 1, there are two directions that 
may align with the force, and, moreover, align- 
ment may not occur at all. Clearly there is much 
more work to be done to achieve complete align- 
ment in this case, and it is not as clear how to go 
about doing it. Furthermore, we argued at the 
end of the introduction that a non-aligning ob- 



ject may be modified so that it becomes axially 
aligning. 

We also make an assumption about the rota- 
tional symmetries that T possesses. Actually, 
the assumption that T be axially aligning al- 
ready puts restrictions on the symmetries the 
object can have. Any symmetry axis of T must 
be an eigenvector. Thus T may have rotational 
symmetries only about v^, and the most symme- 
try T could possibly have is circular symmetry. 
However, we will assume in what follows that T 
has no rotational symmetries. We will consider 
the effect of symmetry in Sec. VI A 

Without loss of generality we take the real 
eigenvalue of T, A3, to be positive. Now we are 
ready to discuss the two possible forcing pro- 
grams. 

A. Step Function Forcing 

We begin by discussing the simplest possible 
time-dependent forcing program, a mere shift in 
the forcing direction by an angle 9: 

- \xsm9 + zcos9 if i < 0, 

^^*^~ [z ifi>0, ^ ^ 

where 9 — termed the rocking angle — is a fixed 
parameter which gives the angle between the ini- 
tial force and the final force. We have chosen the 
axes such that the force rotates onto the z-axis at 
time zero. In this subsection, we will show that 
this forcing program, if repeated sufficiently with 
an appropriate choice of 9, causes an ensemble 
to become completely aligned. 

Let us briefly describe the argument we will 
present for this claim. We start with an ensem- 
ble axially aligned with a force, and we switch 
the direction of the force. After we do this, the 
ensemble will exhibit a transient response where 
it aligns with the new axis. After this transient 
has passed each object simply rotates around the 
new force at constant angular velocity. Since, 
as stated above, no alignment can occur once 
the steady state motion begins, we are interested 
only in the transient motion. This motion causes 
some orientations to bunch together and others 
to move further apart, but we will show that 



on average the ensemble becomes more ordered. 
Thus after many iterations complete alignment 
is achieved. 

This subsection is organized in the following 
way: first, we introduce a formalism to represent 
the effect of the transient motion on an individ- 
ual object. Next we introduce a way of repre- 
senting an axially aligned ensemble, and we de- 
termine the effect of the transient motion on the 
probability distribution of orientations. Then we 
introduce a way to quantify the disorder in the 
ensemble, and we show the transient motion de- 
creases this disorder on average, and we show 
that this means the ensemble must become or- 
dered. 

Let us discuss the effect of the transient. For 
the sake of concreteness, we will temporarily 
consider the parameter 9 to be fixed. Since the 
transient motion may be complicated, our strat- 
egy will be to think of it as a black box: it ac- 
cepts as input an object aligned with the initial 
force, and outputs an object aligned with the 
final force. Since the dynamics are determinis- 
tic, any two objects starting with a given orien- 
tation must end up in same orientation. Thus 
the black box defines a function from the space 
of initial orientations to the space of final ori- 
entations. Since the objects must initially be 
aligned with F, they have only one orientational 
degree of freedom, and the space of initial orien- 
tations is the unit circle, denoted S^. Similarly 
the space of final orientations is also S^ . To rep- 
resent the input orientation, we first arbitrarily 
pick a member of the ensemble as a reference 
object. Then any object in the ensemble can be 
realized by rotating the reference object by an 
angle (j) about F. We use this angle (p to refer to 
the initial orientation. We can analogously rep- 
resent the output orientation, which we will call 
ip. For consistency, we pick the reference object 
for ip to be the final orientation of the initial ref- 
erence object. With this notation in place, we 
can now represent the effect of the transient by 
the function ip{4>)- We note that our choice of 
initial and final reference objects is summarized 
by the relation ip{0) = 0. If we allow the angle 9 
to be chosen freely again, then the function ^{(p) 
will depend on the parameter 9. We indicate this 
dependence by writing the function as ^jJg. Also, 



it will be convenient to require that the choice 
of initial reference orientation be smooth in 6. 

For our purposes, we could forget everything 
about the object, including T, except for ipg. A 
few general properties of ^g will prove useful. 
First, tpQ is the identity function, since ^ = 
corresponds to a constant force, so each orienta- 
tion advances by the same amount. Because the 
choice ipg{0) = 0, this amount must be a multi- 
ple of 2n. This trivial behavior occurring when 
9 = cannot be used for alignment. 

The simple behavior of ipo{(p) leads to sim- 
ple behavior for small, non-zero 9. First, tpg is 
monotonic for sufficiently small 9. To prove this 
statement, we first give the reason why, in (7), 
we chose the initial force to depend on 9 instead 
of the final force. We made this choice because 
we wanted the combination {9, 4>) to represent an 
initial condition for (6): the object with orienta- 
tion (f> when the rocking angle is 9. This initial 
condition is smooth in {9,(j)), because we chose 
the reference orientation to depend smoothly on 
9. Then since the solution of a differential equa- 
tion depends smoothly on initial condition [15], 
ipe{4') depends smoothly on {9,cj)). Then so does 
the (/)-derivative. Since this derivative is identi- 
cally 1 for = 0, it must be positive everywhere 
for sufficiently small 9. However, for certain ob- 
jects and large enough 9 it may be that Ve is 
non-monotonic . 

Our second result is that as the input orien- 
tation (j) is increased from to 27r, the output 
orientation ip0{(j)) undergoes a net increase of 27r 
as well. This is the proper extension of what we 
found for 9 = 0; it holds for any 9. Mathemati- 
cally, this assertion can be written 



dipe 



2tt, 



(8) 



Since Vo is the identity function, (8) is clear in 
the case 9 = 0. Our strategy to prove (8) in 
the general case 7^ is to continuously con- 
nect this situation with the constant force situ- 
ation. Since the solution of a differential equa- 
tion depends continuously on initial condition, 
the function ipg can be deformed into ipo simply 
by decreasing 9. Therefore these two functions 
must wind the same number of times around S^. 
From the 9 = case, we see that it must wind 



exactly once. This proof of (8) is illustrated in 
Fig 3. This concludes our discussion of how the 
transient acts on individual objects. 

Now we describe how the transient acts on 
the ensemble as a whole. We will characterize an 
axially aligned ensemble by a probability density 
function (pdf) p^fp) which gives the probability 
of a randomly selected object having the orien- 
tation (p. Complete alignment means that this 
pdf is a delta function. 

We start by considering an operation more 
simple than ipg. We will study a uniform shift in 
orientation (henceforth, a "shift"), which can be 
effected by allowing an axially aligned ensemble 
to rotate for some amount of time. Although, as 
we have already stated, a shift does not increase 
alignment, it is be an essential ingredient in the 
forcing program we propose. The orientation 4> 
after the shift depends on the orientation (p be- 
fore the shift as follows: (p = (p + a. The new 
pdf is then given by p{(p) = p{(p)', since the new 
pdf is essentially identical to the old pdf, we see 
that no ordering has been achieved. 

Now let's consider ipg, the effect of a single 
shift in forcing. It is convenient to assume that 
the derivative of ipg is everywhere positive, which 
can be assured, as explained above, by making 
9 sufficiently small. We make this assumption 
for two reasons: it excludes certain pathological 
ipe^s, and it simplifies the transformation law of 
the pdf. To help motivate this assumption, we 
consider a pathological V'ejit is shown in Fig. 4 
(a). Note that while this ipg is a valid example 
insofar as it satisfies the winding number con- 
straint, it is nevertheless clear that no ordering 
can be achieved, because any localized distribu- 
tion will be broadened by the action of ipQ, as 
shown in Fig. 4 (b). The assumption that ipg 
is monotonic, together with the condition that 
Ipg winds only once around S^ , guarantees that 
there is a limit to how much stretching ipg can 
do. 

Now we are ready to examine how the pdf 
transforms under the action of ip (we no longer 
concern ourselves with the value of the parame- 
ter 9 and we write simply ip in place of ipg). In 
fact, it is no more complicated to find how the 
pdf transforms under the composition of ip with 






FIG. 3. Illustration of the proof that tpe wraps once around S^. Four graphical representations of ipe are 
shown for four different values of the rocking angle 0. The color of a point on the curve indicates the initial 
angle 4>, and the azimuthal coordinate of a point on the curve is equal to "ipei^f^)- The radial coordinate has 
no physical significance; its purpose is to avoid self-intersections of the curve. For the far left figure, 6* = 0, 
so that ipe is the identity function. As 9 is increased from 0, the azimuthal coordinate of each point changes, 
deforming the curve. However, since this deformation is smooth, it is impossible to change how many times 
ipg winds around S^ just by changing 6. 
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FIG. 4. (a) Plot of ip{(j>) vs </> for a non-monotonic V' that cannot possibly lead to ordering, (b) A plot of 
a narrowly peaked distribution p((/)), shown in black, together with p(0), the result of acting by ip on p(0), 
shown in red. p is broader and therefore represents a less ordered distribution. 



a shift. The orientation after this operation is 



<P = ^(0 + a), 



(9) 



where a is the the size of the shift; it gives addi- 
tional freedom which may be used to help align- 
ment. From the rule p{(j))d(j) = p{(f))d(f), we find 
that the pdf transforms according to the follow- 
ing formula: 



p{<i>) 



p{4>) 



V''(<^ + a)' 



where (f) = 'ip~^{cp) — a. 

Now that we know how the ensemble can be 
manipulated, we must ask if it is becoming more 



or less ordered. We will use the information- 
theoretic entropy, H [16], to quantify the disor- 
der. Given a pdf, p{4>), the functional H[p] is 
defined as 



H[p] 



p{(j)) lnp{(p) d(j). 



(11) 



(10) 



To justify our use of the entropy, let us recall a 
few of its properties. First, a change in scale in 
say from radians to degrees, causes a change 
in H by an additive constant. In this sense, the 
zero of entropy is arbitrary, similar to the case 
of energy. Another important property is that it 
is maximal for a uniformly distributed variable. 



This means there is an upper bound on the en- 
tropy. However there is no lower bound. In fact 
a pdf hmited to some region of size Ai;^ can have 
an entropy of at most log(Ai?!)). This upper hmit 
makes it clear that the entropy goes to — cxj as 
the distribution becomes more and more narrow. 
These observations lends credence to the inter- 
pretation of i? as a measure of disorder. 

We now calculate the change in entropy re- 
sulting from the transient. The entropy Ha af- 
ter the transient with shift a is given by Hq, = 
— fgip{(l))ln{p{(p)) d(p. Upon substituting (10) 
and changing the integration variable from (f) to 
6, we find H = H + AHn where 



AH^ 



p{cP)ln{iP\cj) + a)) d<p. (12) 



If the new distribution is to be more ordered 
than the old distribution, AH^ must be less than 
zero. However, AHa is certainly greater than 
zero for some choices oip{(j)) and a. For example, 
iip{(p) is strongly concentrated in a region where 
^l^' > 1, then AHq is positive. 

Even though the change in entropy is not 
negative for all values of a, we now show 
that if a is chosen randomly, the expected 
entropy change will be negative. The a- 
averaged change in entropy {AHa) is given by 
^ fgi AHada. Using (12) for AHa, we find 
{AHa) = ^§p{(t))d4>§ln{ij'{u))du. The inte- 
gral over (j) is unity, so any dependence of the 
average entropy change on the initial pdf disap- 
pears. Also, because the log is a concave func- 
tion, ^§ln{tP'{u))du < ln{§ ^^P'{u)du) = 0, 
with equality only in the case that ^ is the iden- 
tity. Thus we conclude (AHa) is negative and 
independent of the initial probability distribu- 
tion p. 

With this information about (AHa) we can 
infer the change in H after many repetitions 
of the rocking procedure with randomly chosen 
shifts. We denote the average change in H after 
many such iterations as AH. The average of AH 
after n iterations with a particular randomly 
chosen sequence of shift angles {ai, a2) •■•j"n} 
is given by 



^h=Iy.^^-'^p^ 



where Pi{4>) is the probability distribution before 
the ith iteration. The ensemble average of AH 
over the (randomly chosen) a^ is then given by 



1 "■ 



(14) 



i=l 



Though Pi in each term of this sum depends on 
previous aj's, it is independent of Oj. Thus each 
term average is given by the {AHa) calculated 
in (12). This average is independent of the pi. 
Thus 



{AHa) = {^Ha 



1 

2^ 



hitp' {u)du (15) 



(13) 



i=l 



Since the expectation of the iteration average is 
the same as the expectation of a single itera- 
tion, H is expected to decrease indefinitely after 
many iterations. Moreover, since tp' has an up- 
per bound and is bounded below by 0, the range 
of values AHa can take at each iteration is also 
bounded. Then the central limit theorem [17] 
guarantees that AH approaches {AHa) as the 
number of iterations n goes to infinity. 

Now we see why it is important to compose 
tjj with the shift. With no phase shift, it is dif- 
ficult to rule out the possibility that upon be- 
coming sufficiently ordered, the system would 
consistently lose its order at the next iteration. 
However, with randomness it is impossible for 
the system to consistently pick shift angles that 
cause disorder. 

The indefinitely small entropy found above 
does not guarantee alignment to a single orien- 
tation: there may be two (or more) final orienta- 
tions. However, a recent argument [18] suggests 
that the final state has only one. This concludes 
our proof that repeatedly applying (7) with suffi- 
ciently small 9 will cause the ensemble to become 
completely aligned. 



B. Rotating Force 

One shortcoming of the approach described 
in the previous section is that the final orien- 
tation is arbitrary. In this section we consider 
an alternative forcing program, in which the fi- 
nal orientation is controlled. In this program, 
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the force rotates about the z-axis with constant 
angular velocity: 

F{t) = z cos 00 + {x cosilt+y sini^t) sm6Q. (16) 

Here, 9 and il. are adjustable forcing parameters. 
Without loss of generality, we take fi > and 
< 9 < TT. Since we have not specified the 
direction of the z-axis, it is no more general to 
allow the force to rotate about an arbitrary axis. 
This axis will be indicated by the vector $7, which 
points in the direction of rotation of the force 
and has length 0. 

If the ensemble becomes completely aligned, 
all objects will undergo an identical motion. 
Since the force is not constant, this motion could 
be very complicated. However, we will show in 
this section that for a certain choice of parame- 
ters, the ensemble does become aligned, and the 
subsequent motion is simple. This subsection is 
organized as follows. First, we will describe the 
simple motion occurring after alignment. Next, 
we will examine what choice of forcing param- 
eters gives rise to this simple motion. Finally 
we will show that with this choice of forcing pa- 
rameters the motion is locally stable. Numerical 
studies reported in Section V A indicate that this 
motion is globally stable. 

We first describe the simple motion which will 
occur after alignment. This motion is rotation 
with the same constant angular velocity as F: Q. 
The direction of the force relative to the object 
does not change, so the state of the system at a 
later time is just a rotated version of the initial 
state, and the motion persists. Since the object 
and the force are rotating together, we call this 
motion "co-rotation". It is characterized by the 
condition il = TF, called the co-rotation condi- 
tion. Notice, though, that this motion may not 
be possible for all choices of parameters. For ex- 
ample, cannot be larger than the maximum of 
||TF||. 

We now determine explicitly what conditions 
on the forcing parameters are necessary for there 
to be a co-rotating state. Here we shall express 
F, T, etc. in the body frame, indicated by a sub- 
script "b". We first observe that for any initial 
Fh, there is a unique O;, that gives co-rotation, 
namely fif, = tJf, = T^Ff,. Next, the length of this 



Qh and the angle it makes with the initial Fb de- 
termines forcing parameters ft and 9 compatible 
with this Fh- Thus each F}, generates a particu- 
lar choice of the parameters fi and 9. However, 
in practice, it is not the direction of the force, 
but the parameters O and 9 which are directly 
controlled. For a given and 9, there may sev- 
eral FbS which generate these parameters, i.e., 
several co-rotating orientations of the object in 
the lab frame. On the other hand, it may be the 
case that there are no such .Ffe's. We will find 
a particular choice of and 9 such that there 
is certainly at least one Ff, which generates the 
chosen fi and 9. 

We begin with an object with twist matrix 
T(,, and apply a force Fh- In this frame, the co- 
rotation condition is ThF^ = 0,^. So 0,^ is fixed 
by our choice of Ff,. The parameter Q is fixed by 
the "O constraint" 

Q^ = Qh-^b = F^TlThPh, (17) 

and 9 is fixed by the "^ constraint" 



cos 9 



Fb-ftb 



Fwm 

llTbFbll 



(18) 



Now that we know what forcing parameters 
are generated by a particular choice of F^, we 
will ask what range of forcing parameters could 
be generated by some F^. First, let's consider 
the il constraint, (17). Since T^T^ is a sym- 
inetric, positive definite matrix, its eigenvalues 
are real and positive. Therefore, the range val- 
ues that FffTjTi,Fi, takes is the interval [S3, Si], 
where S3 (Si) is the smallest (largest) eigen- 
value of T^Tf). Evidently Q,"^ must lie in the 
range S3 < f]^ < Si for co-rotation to be pos- 
sible. If some value of Q is chosen, then Fj, is 
constrained to lie on the corresponding level set 
of F^T^T^Fi). A level set is the intersection of 
the ellipsoid defined by (17), and the unit sphere 
since we took all forces to have unit length. We 
will call each one of these level sets "0-curves" . 
Notice that the Q-curves are symmetric under 
reflection through the origin, and they are gener- 
ically a pair of closed loops (hence the plural) . 

Let's examine how the O-curves change as O^ 
changes. If 0^ = S3, then the only F^'s capa- 
ble of producing such a small Q point along the 
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eigendirection associated with S3, so that the 
i7-curves are simply a pair of antipodal points. 
As ri^ increases from S3, the f]-curves form in- 
creasingly larger loops around (±) the eigenvec- 
tor corresponding to S3. As il^ reaches the mid- 
dle eigenvalue of T^T^, the two loops meet at 
(lb) the corresponding eigenvector, and as Q,"^ in- 
creases more, the loops shrink to (it) the eigen- 
vector corresponding to Si. The il-curves look 
similar to the orbits of the angular momentum 
in the body frame of a freely rotating object as 
discussed previously. 

Let's proceed analogously with the 9 con- 
straint (18). For some choice of 9, we call the set 
of Ffe's that solve this equation ^-curves. First 
we consider the possible range of 9. We know we 
can achieve 9 = Ohy choosing F^, = v^. However, 
the other extreme, ^ = tt is unattainable since T 
has no negative eigenvalue. So the range of val- 
ues 9 can take goes from zero up to some maxi- 
mum, which is less than vr. Notice that, like the 
il-curves, the ^-curves are symmetric with re- 
spect to reflections through the origin. Let's see 
how the 0-curves change as 9 is changed. When 
= 0, the ^-curves are just the pair of points 
±1/3. As 9 is increased, the ^-curves expand from 
these two points. Below we shall consider only 
the regime of small 9. 

For a choice of forcing parameters to be com- 
patible with co-rotation, there must be an F^ 
that simultaneously meets the conditions (17) 
and (18) — that is to say the r2-curves and the 9- 
curves must intersect. Phrasing the problem in 
terms of the intersection between Q, curves and 9 
curves leads to an important observation. If one 
of the curves intersects one of the 9 curves, 
then, generically, they must intersect again. Ad- 
ditionally, by inversion symmetry, the antipodes 
of these to intersection points will also be in- 
tersection points. Therefore, a choice of forcing 
parameters compatible with co-rotation actually 
gives rise to four co-rotating states. An example 
of r^-curves and ^-curves is illustrated in Fig. 5. 
This concludes our discussion of which forcing 
parameters give rise to co-rotation. 

The next topic to be discussed is the stability 
of the co-rotating state. Each of the four co- 
rotating states may have a different stability. For 
there to be a globally stable co-rotating state, 




FIG. 5. The Vl- and 6'-curves. We have chosen an ar- 
bitrary twist matrix, and chosen Vi and 9 to be com- 
patible with co-rotation. This figure shows a part of 
each rj-curve (dashed) and one of the 6'-curves (sohd). 
For reference, v^ is shown by a black dot. Sections 
of the two O curves are visible above and below the 
dot. Also, there is another 0-curve antipodal to the 
one visible in the figure. We see that the there are two 
visible intersection points of these sets of curves, and 
by the inversion symmetry, the antipodes of these two 
intersection points are also intersection points. Thus 
there are four co-rotating states. 



one and only one of these states must be stable. 
We will show this to be true. 

To this end, we will analyze the object's mo- 
tion in the frame rotating with F at constant 
angular velocity Q. Note that in general, the 
object need not co-rotate with F and thus the 
object frame is in general distinct from the rotat- 
ing frame. We define T/ to be the twist matrix in 
this frame, so that T/ is defined by the following 
formula: 



T 



W/ 



(19) 



where M(t) = exp[r2^t]. Notice that if the object 
is simply rotating with angular velocity fi, then 
T/ is a constant. We now find the equation of 
motion for T/. Applying (19) to (6), we get, 
after some algebra. 



t/ = [(T,Fo-!^)^T/], 



(20) 



where Fq = F{Q). This equation has the same 
form as the lab frame equation of motion (6). 
The only difference between these equations is 
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that the angular velocity TF has been replaced 
by an apparent angular velocity TjFq — Q, as 
is appropriate for a frame rotating with angular 
velocity 0. Studying (20) reveals that there is a 
fixed point if the co-rotation condition, TjFq = 
Q, is satisfied. This was the purpose of choosing 
the rotating frame. 

We now study the stability of these fixed 
points. To this end, we fix Fq and il. and as- 
sume that the twist matrix Tj is a fixed point 
(i.e. TfFo = fi). We will perform a linear stabil- 
ity analysis of this fixed point. We may parame- 
terize the departure of the twist matrix from Tj 
using an infinitessimal rotation vector fj. We put 



Ti = Tf + [ff'',Tf]. 



(21) 



Next we obtain the linearized equation of mo- 
tion for fj. We use the form (21) in (20) keep- 
ing terms linear in fJ. In the right hand side 
of (20), the first argument of the commutator 
represents the apparent angular velocity in the 
rotating frame; it must be first order in rj. Since 
we are only working to first order, we may take 
only the zeroth order part of the right argument 
of the commutator. We find 



[^^^T/] 



■■,Tf]For,Tf]. (22) 



Eq. (22) shows that {ff - [77^,T/]Fo)'' com- 
mutes with Tf. Unless the former is zero, this 
means that Ty is unchanged upon an infinitesi- 
mal rotation in the ff — [f]^,Tf]Fo direction. In 
this case, Tj has an axis of symmetry. Since we 
have restricted our attention to asymmetric T's, 
there can be no such axis; thus ff must equal 

Further manipulation yields 



= 7/xT/Fo-Tr/(7/xFo) 
= -TfFoxf] + Tf{FoxfJ) 

= -iTfForff+TfF,^ff, 



(23) 



and so we have the following differential equation 
for fj: 



fJ=[TfF,^-{TfFor]ff. 



(24) 



From this equation we find that the condition for 
the fixed point Tj to be locally stable is that all 



the eigenvalues of the "stability matrix" TjFq — 
(TfFo)^ have negative real part. 

We have seen above that a given rotating 
force gives in general four possible co-rotating 
orientations of the object. Our goal of complete 
alignment would require that only one of these 
orientations survives after a long times. While 
this does not happen for all values of forcing pa- 
rameters, we will now show that there is a suit- 
able choice of forcing parameters for which a sin- 
gle stable fixed point exists. 

We consider the regime where 6 is very close 
to and $7 = A3. Let's determine where the fi- 
and 9- curves intersect in this case. Since 9 is so 
small, the 0-curves are very small circles around 
zbiTs. The curves are a little more complicated. 
Since we have chosen 0, to be the one generated 
by Fh = fj-s, fJs is on the fi-curve. Moreover, one 
readily verifies that the Jl-curve must be smooth 
at U3 as follows. One can show from the form 
(25) of T given below that for the asymmetric 
objects considered here, V3 cannot be an eigen- 
vector of T"^, and thus it cannot be an eigenvec- 
tor of T-^T. Since the i7-curves are degenerate 
only at the eigenvectors of T T, the fi-curves 
must be smooth at t% as claimed. Thus, at fJ^ 
the Q-curve is almost straight on the scale of the 
9-curve. Therefore, there will be exactly four in- 
tersection points, two near v^ and two near —V3. 
This situation is shown in Fig. (6). 

It may be illuminating to note the similarity 
between the force used in the forcing program 
proposed here, and the magnetic field used in nu- 
clear magnetic resonance (NMR) . In NMR, there 
is a large DC magnetic field applied, causing the 
spins to have a response characterized by some 
angular velocity cJ. In the plane perpendicular 
to the DC magnetic field, a small (corresponding 
to ~ 0) AC magnetic field is also applied with 
angular frequency fi = w. 

Now, we will show that if the forcing param- 
eters are chosen in this regime, only one fixed 
point is locally stable. More specifically, we will 
show that there is only one fixed point at which 
all the eigenvalues of the stability matrix have 
negative real part. To do this we will set = A3, 
and pick for the parameter 9 some arbitrary, 
small value 9*. We fix our attention on just 
the four fixed points resulting from this choice 
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FIG. 6. Sphere of Fb's for the case il = A3 and small 
9. As in Fig. 5, the dashed curve is the i7-curve, the 
solid curve is the 0-curve, and the black dot indicates 
the position of v^. The fJ-curve passes through U3, 
since fi = A3, and the 0-curve is a small loop around 
■1)3, since was chosen small. As seen in the fig- 
ure, this means that these curves intersect in exactly 
two places. (The other two intersection points are 
the antipodes of the ones shown in the figure.) The 
light region is the "stable region" , where co-rotating 
states are locally stable. Only one of the visible in- 
tersections lies in the stable region; both of the fixed 
points that can't be seen are outside of the stable 
region, as explained in the text. 



of forcing parameters. 

First we will consider the two fixed points 
where F^ is near —v^. Since these F^s are so 
close to — ^3, we would expect the behavior to be 
similar to the case where F^ is —V3. Therefore, 
let us study this latter case. Since —V3 is an 
eigenvector, Q is parallel to F^, so 9 = and 
the force is constant. We know from previous 
work [1] that, at this fixed point, two directions 
are unstable while the third is neutrally stable. 
The stability matrix varies smoothly with 9 near 
9 = 0. Thus for a sufficiently small value 9* of 
9, the unstable eigenvalues remain unstable, so 
that these fixed points are unstable. 

We may perform a similar analysis for the 
two fixed points where Ff, is near +V3. Here, 
for ^ = 0, there are two stable eigenvalues and 
a zero eigenvalue. Therefore, the two 9 = 9* 
fixed points must be stable in two directions, but 
overall stability requires that they also be stable 
in the third direction. 



To determine the stability in the third direc- 
tion, we will consider the 9 = 9* fixed points as 
perturbations of the ^ = fixed point. When 9 is 
increased from 0, the co-rotating Ff,'s move out- 
ward from ^3 along the O-curve in opposite di- 
rections. The stability matrix in the body frame 
is TfoFj^^ — (TftFft)^; it is linear in F^. Since the 
change in Fi, is opposite for the two fixed points, 
the changes in the stability matrix, and, there- 
fore, the eigenvalues are opposite. The neutral 
eigenvalue either has opposite signs at the two 
fixed points or it remains zero for both. All that 
is left to do then is demonstrate that the change 
in the eigenvalue is non-zero. 

Our strategy is to assume this change is zero 
and obtain a contradiction. First, we will choose 
a coordinate system, then we will write out the 
eigenvalue equation to lowest order in 9. We be- 
gin by picking our body frame coordinate axes, 
X, y, and z, so that 



T. 



hi ti2 


' 


^21 *22 





b 


A3 



(25) 



Then the tangent direction of the 0-curve at 
^3 is the direction in the tangent space of the 
unit sphere which is perpendicular to TJTi)V^ = 
(0,6A3,A|). This direction is x. Working to 
first order in 9, we write Fjf = (9,0,1). The 
third eigenvector of the stability matrix will have 
changed by some small amount w = {wi,W2,0) 
so that the eigenvector is now z + w, where 
wi and W2 are small when 9 is small. Since 
we assume the eigenvalue is zero to lowest or- 
der, the eigenvalue equation satisfied by w reads: 
(Tfozx - {Thz)'')w + {Tb9x'< - {Th9x)>')z = 0, 
where we have dropped terms higher than first 
order in 9. Notice that this equation is linear in 
9, wi, and W2, and so we may, by standard meth- 
ods, determine if this equation has non-trivial 
solutions. The result is that the equation has 
non-trivial solutions only when at least one of 
the two conditions are satisfied: (1) tn = A3 and 
^21 = or (2) 6 = 0. The case (1) implies that A3 
is an eigenvalue of T with multiplicity 2, contra- 
dicting the fact that T has complex eigenvalues. 
The case (2) implies that the object twist matrix 
has two-fold rotational symmetry about the z- 
axis, which we have assumed not to be the case. 
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We are forced to conclude that there is no soki- 
tion to the eigenvalue equation if the eigenvalue 
is taken to be zero. Evidently, then, the eigen- 
value cannot remain zero to linear order. This 
implies, as argued above, that there is only one 
stable fixed point. This concludes our proof that 
there is one locally stable fixed point if the forc- 
ing parameters are chosen so that ^2 = A3 and 9 
is small. In Section VB, we will discuss numer- 
ical simulations showing that the fixed point is 
in general globally stable. 



be discussed in more detail for each forcing pro- 
gram. (E) The orientation of each matrix about 
the z-axis is determined. (In the case of the ro- 
tating force, a constant force along the z-axis is 
first applied to ensure that the twist matrices 
are at least axially aligned.) (F) This process is 
repeated for many different random initial twist 
matrices in order to explore the range of different 
possible outcomes. In the following two subsec- 
tions we will discuss step (D) of the simulations 
and present results both forcing programs. 



V. NUMERICAL SIMULATIONS 



A. Step Function Forcing 



We carried out numerical simulations to test 
the results of applying the forcing programs out- 
lined in Sections IV A and IV B above. In addi- 
tion to confirming the behavior predicted in the 
previous section, the goals of the simulations are 
to get a sense of: (1) whether the decrease in en- 
tropy resulting from the step function forcing is 
indeed accompanied by complete alignment; (2) 
whether the condition of having ip monotonic is a 
necessary condition for achieving complete align- 
ment; (3) whether the locally stable fixed point 
of the rotating force program is globally stable; 
and (4) whether alignment is possible outside of 
the small 9 regime. The results of our simulation 
will be presented in this section and the above 
four points will be discussed. 

The simulations were done the following way: 
(A) we obtain a twist matrix, T, by generating a 
3x3 matrix with entries drawn uniformly from 
the interval [—1, 1]; we reject any symmetric or 
non-aligning matrices to ensure that the matrix 
we obtain is axially aligning. (B) We subject 
the twist matrix to a constant force along the 
z-axis so that its real eigenvector is aligned with 
that direction. This is done by solving the dif- 
ferential equation (6) with T as initial data. (C) 
500 angles are drawn uniformly from the interval 
[0,27r], and new twist matrices are then formed 
by rotating the twist matrix found in step (B) by 
the random angles. These 500 T's form our axi- 
ally aligned ensemble. (D) These new twist ma- 
trices are then subjected to a forcing program, 
which again is done by numerically solving (6) 
with the appropriate initial conditions. This will 



As mentioned in the previous paragraph, our 
aim is to test if complete alignment occurs, since 
this is not guaranteed by the decreasing of en- 
tropy. We also test whether it is possible to ef- 
fect complete alignment even when V'(</') is non- 
monotonic. To study the behavior, we use the 
program described in the previous paragraph, 
where step (D) is carried out in the following 
three stages: first the objects, formed by step 
(C), which are aligned with the z-axis, are sub- 
jected to a constant force along the rc-axis until 
they are aligned to a resolution of 10 milliradi- 
ans. This stage is the physical manifestation of 
applying ip. Second, they are all allowed to ro- 
tate for an additional time randomly drawn from 
a uniform distribution ranging from zero to the 
period of rotation. This stage applies the shift. 
Finally an explicit 90° rotation about the y-axis 
is applied to each T matrix, so that they are 
again aligned with the z-axis. We then repeat 
the process starting from stage 1. Each repeti- 
tion is a "step function iteration." The step func- 
tion iteration is illustrated in Fig. 7. We subject 
the system to as many step function iterations as 
desired (which in our simulations was 100). Af- 
ter each iteration, the orientation of each object 
is determined, and the entropy of the ensemble 
is found using a nearest neighbor estimate [19]. 

Some of our results are shown in Fig. 8. Each 
row of plots corresponds to a different twist ma- 
trix. The first row is associated with a twist 
matrix whose ip function is monotonic. From 
left to right, the plots in this (and every other) 
row are: ip{4>) vs (p, the final orientation after 
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FIG. 7. Above is a schematic illustration of one step function iteration. On the far left of the figure is the 
lab coordinate axes. There are three stages shown, and for each stage we show the force vector in green, 
and a representative object of the ensemble in red (in reality there are many objects in the ensemble, but 
only one is depicted for ease of viewing). At the beginning of an iteration, the ensemble is axially aligned 
with the z-axis, and a force is applied along the z-axis. After the ensemble is given some fixed time to 
evolve, the ensemble becomes axially aligned with the a;-axis. After this, the ensemble is allowed to rotate 
for a random amount of time. The next step in the simulation is to rotate the objects by hand so that they 
become aligned with the z-axis again. After this step, we are ready to begin the next iteration. 



all 100 iterations vs initial orientation, and the 
entropy (in bits) as a function of iteration num- 
ber. We can see from the second plot in the first 
row that the final orientations all share a com- 
mon value; i.e., the objects have aligned. Thus 
we did not run into the problem of having mul- 
tiple final orientations, even though this would 
have been consistent with our proof of indefi- 
nite decrease in entropy. In the third plot of the 
first row, there is a roughly linear trend between 
entropy and iteration number, consistent with 
our prediction in (12). The final entropy, which 
is approximately —16 bits, can be attributed to 
numerical error. In our simulations, each object 
that had a monotonic -0 exhibited behavior simi- 
lar to this example. In particular, we found that 
complete alignment was attained. 



The next row shows that ordering can occur 
even with a non- monotonic "0- We see in the 
second plot that although ip is non-monotonic, 
all the objects have the same final orientation. 
The entropy decrease is more erratic in this case. 
However, not all twist matrices that we exam- 
ined gave alignment. For example, the twist ma- 
trix of the third row resulted in a random dis- 
tribution of final orientations. Correspondingly, 
we see that the entropy remains near its initial 
value. 



B. Rotating Force 

In the case of the rotating force, our goals 
were to test if the stable fixed point is always 
globally stable when 9 is small, and to deter- 
mine the behavior for larger 9, where we have 
not even shown a stable fixed points exists. To 
test global stability for small 9, we performed 
step (D) of the simulation by starting with an 
axially aligned ensemble, and applying a sedi- 
menting force of the form (16) to an ensemble 
of objects, with 9 small, and ri = A3. To mea- 
sure the degree of alignment of the ensemble at 
any point during the forcing, we would make a 
copy of the ensemble, and apply a constant force 
along the z-axis so the ensemble would at least 
be axially aligned. Then the orientation of each 
object is determined, and the standard deviation 
of the set of orientations is computed. Then rela- 
tive orientations are determined and the ensem- 
ble is considered aligned if this standard devia- 
tion is less than 0.5 milliradians. Having done 
this simulation on 1000 different randomly gen- 
erated twist matrices, we found that alignment 
is always achieved if 9 is chosen to be sufficiently 
small and the system is evolved for a sufficiently 
long time. The typical value of 9 was about 0.07 
radians, but, as we will show, 9 could have been 
chosen to be larger in many cases. The typical 
amount of time required for alignment is about 
500 periods of the forcing. 

However, it was occasionally the case that 
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FIG. 8. Results of simulating the alternating forcing program for three different twist matrices. Each row 
of plots corresponds to a different twist matrix. In each row, the first plot shows V'(0) vs 0; the second plot 
shows (/>finai7 which is the orientation resulting from 100 step function iterations, vs (/)initiai; and the third 
plot shows the entropy as a function of iteration number. The line in the first row indicates the decrease of 
entropy in bits predicted by (12) 



alignment took a very long time. The twist 
matrices which took longest to align were ones 
where v^ was very near one of the eigenvectors, 
call it w, of T-^T. Before determining why this 
leads to slow alignment, we identify experimen- 
tal behavior connected with nearly coinciding 
eigenvectors. In the coordinate system of (25), 
T'^Ws = A3(0, 6, A3). Thus we find that v^ is an 
eigenvector of T T only if either A3 = or 6 = 0. 
Therefore, if W3 is nearly and eigenvector of T'^T, 
either A3 or b must be small. The meaning of a 
small A3 is that the object has a slow response 
to a constant force, and thus the rotating force 
is required to rotate very slowly. The meaning 
of 6 = 0, as can be verified by direct computa- 
tion, is that z is a two-fold symmetry axis of T. 
Therefore, the problem twist matrices are ones 
where either the force is required to rotate very 
slowly or the twist matrix was nearly two-fold 
symmetric. 

The small difference between V3 and w leads 
to slow alignment for two reasons to be explained 
below. First, 6 must be small in order to satisfy 



both the 0, and the 6 constraints. Second, a 
small 9 entails slow alignment. 

We first explain why 6 must be small when 
the eigenvectors nearly coincide. On the one 
hand, the T-^T close to vs may correspond to an 
extremal eigenvalue of T T. In this case, one re- 
curve is closely localized around w. Thus if 6 is 
chosen too large, the 0-curve will completely cir- 
cumscribe the f] curve, and these curves will not 
intersect. Thus 9 must be small if V3 is close to 
w in this case. On the other hand, w may corre- 
spond to the intermediate eigenvalue. Then the 
two solution curves of (17) are almost touching 
at v^. Therefore, if 6 is chosen too large, the 9 
curve will intersect both Q curves, and there will 
be a total of 8 fixed points, possibly two of them 
stable. Thus 9 must be small in this case as well. 

The smallness of 9 slows the alignment be- 
cause as 9 goes to 0, the third eigenvalue of 
the stability matrix also must go continuously to 
zero. Thus the time required for alignment be- 
comes very long, and the numerical simulation 
is difficult. Of course this problem is relevant to 
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physical applications as well as numerical simu- 
lations. 

Next we consider the behavior outside of the 
small 9 regime. To explore the full range of 9 
by simulation, multiple copies of the initial en- 
semble were created, each one subjected to a dif- 
ferent forcing program. The ^'s of these forcing 
programs were evenly distributed between and 
vr, but as before we set i7 = A3. 

We characterize how "well" an object aligns 
as a function of 9 in the following sense. If the 
ensemble is sufficiently converged, then the sys- 
tem is well described by a linear approximation, 
and the standard deviation of the orientations 
must decay exponentially with time. We will use 
this decay rate to quantify how well the object 
aligns. We know that this decay rate is small 
when 9 is small, but we also have reason to be- 
lieve that alignment is impossible if 9 is too large. 
Thus there is some intermediate optimal 0, and 
studying the alignment rate vs. 9 should give an 
indication of what the optimal 9 tends to be. 

Fig 9 shows simulation results for 10 differ- 
ent twist matrices. We parameterize each twist 
matrix by the angle C, between v-^ and the closest 
eigenvector of T-'^T. We argued above that if the 
C, angle for a particular T is small, then 9 must 
be small in order for alignment to occur, and 
the alignment must be slow. Although the align- 
ment for small C, is typically slow, we find that 
for any C, — small or large — the optimal 9 is typ- 
ically larger than (", showing that not only is our 
method robust to picking a large value for 0, but 
choosing a large 9 is often beneficial. Remark- 
ably, the optimal alignment rate can be many 
multiples of the rotational frequency, although 
the total time required for alignment is even then 
typically a few rotational periods, since align- 
ment in the non-linear regime is slower. Addi- 
tionally, we found that when alignment did oc- 
cur, it was always accompanied by co-rotation. 



VI. DISCUSSION 
A. Assumptions 

In the above derivations, we have made many 
assumptions about the nature of the medium. 



the nature of the forcing, and the nature of the 
sedimenting object. We will now take inventory 
of these assumptions and discuss their validity. 
One of the fundamental assumptions we made 
in our approach to this problem was to neglect 
interactions between the different objects in the 
ensemble, thereby showing that the orientational 
alignment demonstrated here does not require 
these interactions. Thus the underlying mecha- 
nism is different from the one encountered with 
weakly coupled dynamical systems such as a row 
of similar mechanical clocks mounted on a wall 
[20]. However, in practice there is a coupling: a 
given object in the suspension experiences forces 
and torques owing to the flow emanating from 
nearby sedimenting objects. These interactions 
may enhance or diminish the phase locking ef- 
fects shown above. In any case, the interactions 
can be reduced by diluting the suspension. 

Being acted upon by a force, the object will 
move. We assumed, as would be appropriate for 
objects of colloidal length scales, that the result- 
ing motion is characterized by a small Reynolds 
number. Typical objects up to the size of a mil- 
limeter easily satisfy this criterion. Experiments 
[1] up to Reynolds numbers of several hundred 
showed the same sedimentation behavior as for 
Reynolds numbers below 1. However, inertial 
effects at higher Reynolds numbers may well al- 
ter our results. The assumptions listed so far 
allowed us to state that the angular velocity of 
an object at any point in time depends only on 
the force applied to it at that moment, and it 
depends on the force linearly through the twist 
matrix. 

In writing (4), we have regarded the twist ma- 
trix in the body frame as constant. Since this 
matrix depends only on the shape and mass (or 
charge if the forcing is electrophoretic) distribu- 
tion of the object, we have essentially assumed 
that the shape of the object does not change 
with time, or, in other words, the object is rigid. 
This assumption is well satisfied for many ob- 
jects. Still, any given object may be deformed by 
the hydrodynamic stresses. A sufficiently large 
deformation creates corrections to the linear re- 
sponse we have treated here. 

We have also assumed that the objects in 
the ensemble have identical twist matrices. In 
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FIG. 9. Rates of full alignment of a sample of ten randomly generated twist matrices subjected to rotating 
forcing programs of various 6 as described in the text. For clarity two views of the three-dimensional plot 
are shown, and color is used to distinguish points corresponding to different twist matrices. The vertical 
axis shows the rate of alignment t~^, expressed in units of the respective object's natural rotation period 
under constant force. The back-to-front axis represents the angle C between the stable eigenvector v^ of T 
and the nearest eigenvector of T-^T. The left-to-right axis is the angle 9 normalized to this (. Left-to-right 
curves show the effect of increasing the angle 9 between the force and the rotation axis for a given twist 
matrix. The interpolating curves are guides for the eye. Points corresponding to values of 9 that did not 
lead to alignment are omitted from the plot. 



practice any small differences in the twist ma- 
trices leads to misalignment. In the case of 
the step function forcing, one might worry that 
small changes is the twist matrix might lead to 
a vastly different aligned state after many itera- 
tions. However, the ability of our simulations 
to work despite representing the twist matrix 
to only machine precision indicates that this is 
not a problem. In the case of a rotating force, 
whether or not the i7-curves and 0-curves inter- 
sect should be insensitive to small changes in 
the object shape. Therefore, we would still ex- 
pect an ensemble of slightly different objects to 
co-rotate with each other, given appropriately 
chosen forcing parameters. 

We have supposed that the object rotates de- 
terministically in response to F. In practice, 
thermal fluctuations add random rotations to 
this response. Our approximation is valid pro- 
vided that the object size is sufficiently large, 
since the size of the thermal fluctuations de- 
creases with increasing object size. The min- 
imum object size will be estimated in Section 
VI B. 



twist matrix was "axially aligning", with only 
one real eigenvalue. The motivation for this as- 
sumption was that under a constant force, these 
object become axially aligned. If this assump- 
tion is not in force, then there are two differ- 
ent axes in the body that may align with a con- 
stant force, so there is no way to characterize the 
transient by a single function ip and the formal- 
ism used to demonstrate alignment completely 
breaks down. However, the rotating force may 
still align, and this is a possible direction of fu- 
ture research. 

A further assumption made above was to con- 
sider generic objects with no rotational symme- 
try. Since T is axially aligning, it necessarily 
has a preferred axis which must be preserved by 
any symmetry. Therefore, the only symmetries 
T can possess are discrete n-fold rotations about 
this axis or complete circular symmetry about 
the axis. With complete rotational symmetry, 
further alignment is impossible since the axially 
aligned T's are already identical. The only n- 
fold symmetry to consider is 2-fold symmetry 
since higher symmetries are impossible. ^ In the 



Another assumption we made was that the ^ A rotational symmetry of a tensor must simultane- 
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case of a two-fold symmetry, the step function 
forcing results are easily extended: where previ- 
ously an objects orientation was represented by 
a number (p in the interval from to 27r, here we 
need only consider orientations in the interval 
from to vr. The function ip maps this interval 
into itself. Since there is no essential difference 
between this interval and the full interval from 
to 2-K, the proof of alignment goes through un- 
modified. However the case of rotating force is 
not so straightforward. Here V3 is necessarily an 
eigenvector of T T, and we are unable to use 
the same arguments that we did in the case of a 
non-symmetric T. This is a direction for future 
work. 

In addition to these global assumptions we 
made further restrictions particular to either 
forcing method. For the step- function method, 
we assumed that the function ■0 was monotonic. 
We found this could be easily attained by mak- 
ing the rocking angle 9 small enough. In the case 
of the rotating force, we found that a good choice 
of parameters for the rotating force program was 
similar to what is done in nuclear magnetic res- 
onance: we have the DC component of the force 
be much larger that the AC component, and we 
have the AC component rotate at the angular 
velocity given by the response of the object to 
the DC component. Although making 9 small 
should be easy to do experimentally, it may not 
be easy to set Q = A3, since A3 may not even be 
known. One way to surmount this problem is to 
gradually increase the frequency with time. As 
the frequency is swept through A3, the ensemble 



ously be a symmetry of both the symmetric and an- 
tisymmetric parts. The antisymmetric part has the 
form u^ for some u. The only symmetries the anti- 
symmetric part has, then, is rotations about u. The 
symmetric part defines three eigenvalues Ai, A2, and 
A3, with eigenvectors vi, V2, and €13. If the eigenvalues 
are distinct, then the eigenvectors are uniquely defined 
up to a sign. Any rotational symmetry in the distinct 
eigenvalues case must then be a vr rotation about one 
of the eigenvectors, and no higher symmetry than two- 
fold rotations is possible. If two eigenvalues are the 
same, then the symmetric part has circular symmetry 
about the non-degenerate eigenvector. In any event, 
rotational symmetries are only possible when u is an 
eigenvector of the symmetric part, but then the only 
symmetries possible are two-fold symmetry and circu- 
lar symmetry. 



will become completely aligned. Once the fre- 
quency gets beyond A3, the objects will evolve 
in a similar way since they are in a similar ori- 
entation, so the ensemble will remain aligned. 



B. Thermal Diffusion 

In Section VIA we noted that if the object 
was not large enough, thermal rotational diffu- 
sion would make alignment impossible. Thus 
there is a minimum object size required for align- 
ment to occur. Our goal in this subsection is 
to estimate this minimum size. To do this, we 
will first estimate the time for an object to lose 
its orientation as a function of size. Then we 
will estimate the time it takes for an ensemble 
to become aligned as a function of size in two 
scenarios. The first scenario is when the physi- 
cal origin of the aligning force is centrifugation, 
and the second scenario is when the force is due 
to electrophoresis. First, we calculate the rate at 
which an object loses its orientation. We will call 
this the decoherence rate F. Denoting the typi- 
cal angle an object has rotated from the aligned 
state as Ac^, the criterion for the ensemble to 
be disordered is that A(j) ~ 1. In order to es- 
timate diffusion, we approximate the object as 
a sphere whose diameter 2R is the span of the 
object [21]. The rotational diffusion constant Dr 
for such a sphere at temperature T is given by 
Dr = UbT /{SiirjE?), where ks is Boltzmann's 
constant. Then over short times t, the mean- 
squared angular displacement Ai;^ is given by 
(A(/>)^ = 2Drt. This formula gives a decoher- 
ence rate F ^ 2D^ = ^I^J^. 

Next we will compare the decoherence rate 
to the rate at which the objects become aligned. 
To do this, we must first find the typical angular 
velocity of an object. First we consider the case 
where a centrifuge is providing the sedimenting 
force, the centrifugal force is F = pi^acV, where 
Pb is the absolute value of the object's buoyant 
density (the difference between the objects den- 
sity and the density of the solvent); ac is the cen- 
trifugal acceleration felt by the object; and V is 
the volume of solvent displaced by the object. It 
is convenient to relate V to the typical radius R 
using a parameter a representing how much of 
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the pervaded volume is occupied by the object: 



V 



Att 



aR . To estimate the angular speed, we 



observe that, by dimensional analysis, a; v/R so 
that u = P Q^fjii ■ The constant of proportional- 
ity /3 specifies the degree of rotation-translation 
coupling and depends on the shape of the object. 
In our prior simulations of physical shapes, we 
found that /3 is typically of order 0.01 [22]. 

The expression for the angular velocity be- 
comes w = ° n''"'' . Now that we have found 
the angular velocity we will estimate the align- 
ment rate. As indicated by the large alignment 
rates in Fig. 9, objects typically align in just 
a few periods of rotation. Therefore we will 
define the alignment time to be ta = — with 
7 = 0.1. Setting the product of the alignment 
time with the decoherence rate equal to unity, we 
obtain the following lower limit on the radius: 

R > \l f.""^ ■ Putting a = 1, r = 300 K, 
Pb = 100 kg/m^, and ac = 10^ m/s^, we find 
R > 0.7 ^m. 

As noted in Sec. II A, our alignment mecha- 
nism depends on a linear coupling between a vec- 
tor perturbation and the rotation of the objects. 
Since asymmetric objects rotate in response to 
an electric field [5], our alignment methods may 
be implemented via electrophoresis. When an 
electric field E is applied, an object will move 
with a velocity v = fiE, where the proportional- 
ity constant /x is called the electrophoretic mo- 
bility. The mobility depends linearly on the zeta 
potential ("• Specifically A* = „ where e is the 
permittivity of the solvent. Assuming, as we did 
in the case of sedimentation, that u = f3v/R 
and ta = — J and again requiring that Fta > 1, 
we get the following lower bound on the radius: 

R > J 4-n-B^nEC • '^° obtain a numerical value for 
R, we assume the strength of the applied elec- 
tric field is lO'^ V/m, the C potential is 15 mV 
(as is necessary to stabilize the colloidal suspen- 
sion [2.3]), and the relative permittivity of the 
solvent is 80, as with water. Assuming the same 
values as before for other constants, we find the 
minimum radius is 2 fim. 

In either case, we have found a lower limit 
for the size of the object, as measured by the 
hydrodynamic radius, of order 1 //m. 



C. Applications 

Polymer science [24], molecular biology [25] 
colloid science [26] and beam lithography [27] 
have made it possible to produce identical asym- 
metric micron-scale objects in large quantities, 
e.g., pollen grains [2S]. These objects often oc- 
cur as dispersions in a liquid. In current practice, 
the objects interact individually and stochasti- 
cally with their surroundings, so that their rela- 
tive orientations are not important. However, 
the oriented samples envisaged here offer the 
prospect of coherent responses to their environ- 
ment. They move together in response to forces, 
or chemical gradients. Thus a group of objects 
in part of the sample could be steered to migrate 
together into a desired region. The oriented ob- 
jects scatter waves in the same way. Thus the 
determination of their structure via x-ray diffrac- 
tion is simplified by reducing the orientational 
randomness. Any anisotropic emissions such as 
fluorescent radiation or Soret currents are also 
enhanced by the common alignment of the ob- 
jects. Thus the possibility of tandem control of 
the orientation of the objects in a sample sug- 
gests new ways for characterizing them and ma- 
nipulating them. 



VII. CONCLUSION 

These possible applications give only a sug- 
gestion of the potential uses of programmed forc- 
ing of asymmetric colloidal objects. Their value 
in practice will depend on the quantitative im- 
pact of degrading effects like thermal fluctua- 
tions and object variability mentioned above. 
Another limitation comes from the twist ma- 
trices encountered in the actual objects to be 
manipulated. We have seen above that these 
matrices can be controlled to a degree, but the 
actual scope for such control remains to be ex- 
plored. Our work to improve this understanding 
is in progress. Conceptually, this work broadens 
our notions of how random objects can be or- 
ganized by programmed external fields. In the 
quantum world, the value of programmed mag- 
netic and electric fields has been well appreciated 
as a means for organizing the disordered quan- 
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turn states of atoms and molecules. This work 
shows that similar means can be used to orga- 
nize colloidal objects. There is every reason to 
anticipate further extensions. 
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